%==========================================================================
%  EXTRACT_TVEL_GRID uses an original .mat file that was interpolated from 
%  a Comsol model file and extract just the temperature and velocity
%  fields.  This was done so that I could recalculate the melting model
%  without having to be logged in to the COMSOL server.
%
%  Files are created for T = 1300-1450C and U = 0.5-10 cm/yr
%
%  Written by Mark Behn, Nov 2014 (in Botswana)
%==========================================================================


% Define spreading rates to loop through
U = [0.5, 1, 2, 3, 4, 6, 8, 10]; %cm/yr
sU = {'U05','U1','U2','U3','U4','U6','U8','U10'};

Tm = [1300,1350,1400,1450];
sTm = {'1300','1350','1400','1450'};

indir = '../../meltmods/WH_DMM/RidgeTYVisc_';

% Interpolation Parameters used when extracting T,U,V from Comsol runs
Pmin = 0.2;
Pmax = 30;
Pinc = 0.05;
rho = 3300;
g = 9.8;

dx = 2e3;
dz = Pinc*1e8/rho/g;
x = -200e3:dx:200e3;
z = -[Pmax:-Pinc:Pmin]*1e8/rho/g;
Pb = Pmax:-Pinc:Pmin;
nx = length(x);
nz = length(z);
[X,Z]=meshgrid(x,z);

soutparams = ' IT Iu Iv Pmin Pmax Pinc rho Pb g dx dz x z nx nz X Z';

for m = 1:length(U)
    for n = 1:length(Tm)
        infile = [indir,sU{m},'_Tm',sTm{n},'_dFdP01_BaG14_interp.mat'];
        outfile = ['grid_',sU{m},'_Tm',sTm{n},'.mat'];
        eval(['load ',infile])
        eval(['save ',outfile,soutparams])
    end
end

